For a multi-site MRI study on MS, we want to understand the effect of a clinical measure, EDSS, on imaging metrics. Generally the distribution of EDSS is not normal, but is bimodal normal, where scores cluster around 2(low severity) and 6(high severity) on a scale of 0 to 10. The goal of this notebook is to first group subjects into low and high severity groups based on EDSS scores. We will then run the simulated data through a Logistic Regression to predict group membership.

In order to get the most power behind the coefficient estimate for EDSS, we need to combine data across the multiple sites. When combining data, we aim to calibrate the scanner metrics before correlating with EDSS, rather than running the model with uncalibrated metrics and an addition protocol variable. Previously we found that the calibration method is generally superior to the protocol method under certain estimated conditions (i.e. subject variance, scanner variance, number of healthy controls at each site). The second goal of this notebook is to test a different calibration method based on transforming the simulated brain volume metrics into percentiles.

We simulate data the same way we did in the previous analyses:

In [1]:

%pylab inline
from utils import *

Populating the interactive namespace from numpy and matplotlib

In [2]:

from IPython.html.widgets import interact, interactive, fixed
from IPython.html import widgets
from IPython.display import clear_output, display, HTML
import sklearn.linear_model as lm
import time

In [76]:

from mpld3 import enable_notebook, disable_notebook
from mpld3 import plugins

In [4]:

params = {'n_hc':50, #number of controls
          'n_ms':[250,250], #number of MS patients at each site
          'age_hc_range':(20,52), # Age range of the controls
          'age_ms_range':[(20,45),(20,45)], # Age range of MS patients at each site
          'd_dist':[.65,.65], # Distribution weighting of EDSS scores for each site
          'alpha_a':-3, # Real coefficient of age
          'alpha_d':[-16,-16], # Real coefficient of EDSS score
          'b0':1644, #Brain Volume Intercept (in ccs)
          'sig_b':[57,57], # Brain variation across subjects (in ccs) for each site
          'sig_s':[16,16], # Variability in brain volume measurement across scanners for each site
          'alpha':[0.7,1.3], # Scaling of each scanner
          'gamma':[20,-20], # Offset of each scanner
          'lo_mean':[2,2], # EDSS distribution parameters at both sites
          'hi_mean':[6,6],  #EDSS distribution parameters at both sites
          'lo_sig': [1.1,1.1], #EDSS distribution parameters at both sites
          'hi_sig': [0.8,0.8], #EDSS distribution parameters at both sites
          }

In [5]:

sim = sim_data(**params)

Here we see the Bimodal-Normal distribution of EDSS at both sites.

The parameters across both sites are identical, so we expect the distribution to be very similar.

In [7]:

enable_notebook()
sns.set_style("white")
_,x=subplots(1)
ppl.remove_chartjunk(x,["top","right"])
x.set_title("EDSS Distribution",size=14)
sns.kdeplot(sim["dis_1"],shade=True,ax=x,label="Site 1");
sns.kdeplot(sim["dis_2"],shade=True,ax=x, label="Site 2");
x.set_xlim([0,10])
x.legend();

Next we will show the multivariate relationship between brain volume, age, and EDSS scores. The graphs show contours along which brain volume is constant. We expect the slopes of the contours to be similar but not identical across sites:

In [71]:

_,x = subplots(1,2,figsize=(12,4))

sns.set_style("darkgrid")
df1 = pd.DataFrame(index=range(params["n_ms"][0]))
df1["BV"] = sim["BV_ms_1"]
df1["age"] = sim["age_ms_1"]
df1["EDSS"] = sim["dis_1"]
df1["Group"] = np.array(["Lo","Hi"])[(sim["llr_1"] > 0).astype(int)]
iax = sns.interactplot("age","EDSS","BV", df1,cmap="coolwarm", filled=True, levels=25,ax=x[0],colorbar=True);
#iax.figure.colorbar(iax)

sns.set_style("darkgrid")
df2 = pd.DataFrame(index=range(params["n_ms"][1]))
df2["BV"] = sim["BV_ms_2"]
df2["age"] = sim["age_ms_2"]
df2["EDSS"] = sim["dis_2"]
df2["Group"] = np.array(["Lo","Hi"])[(sim["llr_2"] > 0).astype(int)]
sns.interactplot("age","EDSS","BV", df2,cmap="coolwarm", filled=True, levels=25,ax=x[1],colorbar=True);
x[1].set_title("Site 2");

x[0].set_title("Site 1");

For both sites, we see that for a given brain volume, as age increases, EDSS decreases. But the slopes do not match perfectly. This is because of the difference in scanner scaling factors.

Assigning Groups: Log-Likelihood Ratio¶

The Log Likelihood Ratio (LLR) is defined as:

\(LLR = log(\frac{Pr_{\theta_L}(E)}{Pr_{\theta_H}(E)})\)

Where \(Pr_{\theta_L}\) is the probability density function(PDF) of a normal distribution centered at an EDSS 2 with a standard deviation of 1.1 (the low group), and \(Pr_{\theta_H}\) is the PDF of a normal distribution centered at EDSS of 6 with a standard deviation of 0.8.

We define the groups as:

Low: \(LLR > 0\)

High: \(LLR < 0\)

Below we plot the distribution of log likelihood ratios at each site:

In [9]:

_,x = subplots(1,2,figsize=(10,5))
[ppl.remove_chartjunk(i,["top","right"]) for i in x]
x[1].set_xlabel("EDSS LLR")
x[1].set_title("Site 2",size=14)

x[0].set_xlabel("EDSS LLR")
x[0].set_title("Site 1",size=14)
ppl.hist(x[0],sim["llr_1"],color=ppl.set1[0],alpha=0.7);
ppl.hist(x[1],sim["llr_2"],color=ppl.set1[1],alpha=0.7);
suptitle("EDSS Lo/Hi Log-Likelihood Ratio", size=16);

Next we will show the plots of brain volume, age and EDSS with the Low and High groups:

In [70]:

_,x=subplots(2,2,figsize=(12,12))
x = x.ravel()
llr = sim["llr_1"] 
ppl.plot(x[0],sim["age_ms_1"][llr],sim["BV_ms_1"][llr],marker='o',linewidth=0,
         markersize=10,color=ppl.set2[0],label="Lo",mew=0.5,mec=ppl.almost_black,alpha=0.9);
ppl.plot(x[0],sim["age_ms_1"][llr==False],sim["BV_ms_1"][llr==False],marker='o',linewidth=0,
         markersize=10,color=ppl.set2[1],label="Hi",mew=0.5,mec=ppl.almost_black,alpha=0.9);
ppl.plot(x[1],sim["dis_1"][llr],sim["BV_ms_1"][llr],marker='o',linewidth=0,
         markersize=10,color=ppl.set2[0],label="Lo",mew=0.5,mec=ppl.almost_black,alpha=0.9);
ppl.plot(x[1],sim["dis_1"][llr==False],sim["BV_ms_1"][llr==False],marker='o',linewidth=0,
         markersize=10,color=ppl.set2[1],label="Hi",mew=0.5,mec=ppl.almost_black,alpha=0.9);
x[0].set_ylabel("Brain Volume",size=14)
x[0].set_xlabel("Age",size=14)
x[1].set_xlabel("EDSS",size=14)


llr = sim["llr_2"] 
ppl.plot(x[2],sim["age_ms_2"][llr],sim["BV_ms_2"][llr],marker='o',linewidth=0,
         markersize=10,color=ppl.set2[0],label="Lo",mew=0.5,mec=ppl.almost_black,alpha=0.9);
ppl.plot(x[2],sim["age_ms_2"][llr==False],sim["BV_ms_2"][llr==False],marker='o',linewidth=0,
         markersize=10,color=ppl.set2[1],label="Hi",mew=0.5,mec=ppl.almost_black,alpha=0.9);
ppl.plot(x[3],sim["dis_2"][llr],sim["BV_ms_2"][llr],marker='o',linewidth=0,
         markersize=10,color=ppl.set2[0],label="Lo",mew=0.5,mec=ppl.almost_black,alpha=0.9);
ppl.plot(x[3],sim["dis_2"][llr==False],sim["BV_ms_2"][llr==False],marker='o',linewidth=0,
         markersize=10,color=ppl.set2[1],label="Hi",mew=0.5,mec=ppl.almost_black,alpha=0.9);
x[2].set_ylabel("Brain Volume",size=14)
x[2].set_xlabel("Age",size=14)
x[3].set_xlabel("EDSS",size=14)

ppl.legend(x[2]);
x[2].set_title("Site 2", size=14)
x[3].set_title("Site 2", size=14)
x[1].set_title("Site 1", size=14)
x[0].set_title("Site 1", size=14)
ppl.legend(x[3]);

ppl.legend(x[0]);
#suptitle("Site 1",size=14)
ppl.legend(x[1]);
suptitle("EDSS Severity Groups",size=16);

Calibration by Percentile¶

For each site, we convert the brain volume metric to percentile by calculating the empirical cumulative distribution function (CDF). The percentile is the CDF evaluated at a given brain volume.

In [13]:

from statsmodels.distributions.empirical_distribution import ECDF

In [14]:

def to_percentile(sim,ref_key,key):
    foo = ECDF(sim[ref_key],side="right")
    #get_p = lambda foo,sim,i : foo.y[np.nonzero(foo.x == sim[key][i])[0]][0]
    #p = [get_p(foo,sim,i) for i in range(len(sim[key]))]
    return foo(sim[key])

In [15]:

sim["pBV_ms_1"] = to_percentile(sim,"BV_hc_1","BV_ms_1")
sim["pBV_ms_2"] = to_percentile(sim,"BV_hc_2","BV_ms_2")

In [69]:

_,x=subplots(1,2,figsize=(12,4))
ppl.remove_chartjunk(x[0],["top","right"])
ppl.remove_chartjunk(x[1],["top","right"])
sns.kdeplot(to_percentile(sim,"BV_hc_1","BV_ms_1"),ax=x[0])
sns.kdeplot(to_percentile(sim,"BV_hc_2","BV_ms_2"),ax=x[0])
sns.kdeplot(sim["BV_ms_1"],ax=x[1])
sns.kdeplot(sim["BV_ms_2"],ax=x[1])
x[0].set_title("Percentile BV MS",size=14);
x[0].set_xlim([0,1])
x[1].set_title("Uncalibrated BV MS",size=14);

In [17]:

age_euro = np.array([24,24,24,38,37,32,28,30,32,32,50])
n_euro = len(age_euro)
BV_euro_real = params['alpha_a']*age_euro + np.random.normal(params['b0'],params['sig_b'][0],n_euro)
sns.set_style("white")
# Lets get scanned twice
#BV_euro_real = np.hstack((BV_euro_real,BV_euro_real))
#age_euro = np.hstack((age_euro,age_euro))
#n_euro = len(age_euro)

sim['BV_euro_1'] = params['alpha'][0]*BV_euro_real+np.random.normal(params['gamma'][0],params['sig_s'][0]*1/sqrt(2),n_euro)
sim['BV_euro_2'] = params['alpha'][1]*BV_euro_real+np.random.normal(params['gamma'][1],params['sig_s'][1]*1/sqrt(2),n_euro)
fig,ax = subplots(1)
sns.kdeplot(sim['BV_euro_1'],shade=True,label="Site 1 Euro",color=ppl.set1[0],ax=ax)
sns.kdeplot(sim['BV_euro_2'],shade=True,label="Site 2 Euro",color=ppl.set1[1],ax=ax)
xlabel("Brain Volume")
title("Brain Volume Distribution of EUROTRIP Subjects",size=14);
ppl.remove_chartjunk(ax,["top","right"])

In [107]:

sns.set_style("darkgrid")
fig,x = subplots(1,2,figsize=(12,4))

"""
ppl.plot(x[0],sim["BV_hc_1"],to_percentile(sim,"BV_hc_1","BV_hc_1"),marker='o',
         linewidth=0,label="HC 1",markersize=15,alpha=0.5,mew=0.8,color=ppl.set2[0])
ppl.plot(x[0],sim["BV_ms_1"],to_percentile(sim,"BV_hc_1","BV_ms_1"),
         marker='o',linewidth=0,label="MS 1",markersize=5,alpha=1,mew=0.5,color=ppl.almost_black)"""

ppl.scatter(x[0],sim["BV_hc_1"],to_percentile(sim,"BV_hc_1","BV_hc_1"),
         marker='o',label="HC 1",alpha=0.5, color=ppl.set2[0],s=200)
ppl.scatter(x[0],sim["BV_ms_1"],to_percentile(sim,"BV_hc_1","BV_ms_1"),
         marker='o',label="MS 1",alpha=1,color=ppl.almost_black,s=25)

Y1 = to_percentile(sim,"BV_hc_1","BV_euro_1")
#ppl.plot(x[0],sim["BV_euro_1"],Y1,
#         marker='o',linewidth=0,label="msHC 1",markersize=10,alpha=1,mew=0.5,color=ppl.set2[1])
scatter1 = ppl.scatter(x[0],sim["BV_euro_1"],Y1,
         marker='o',label="msHC 1",alpha=1,color=ppl.set2[1],s=100)

"""ppl.plot(x[1],sim["BV_hc_2"],to_percentile(sim,"BV_hc_2","BV_hc_2"),
         marker='o',linewidth=0,label="HC 2",markersize=15,alpha=0.5,mew=0.5,color=ppl.set2[0])
ppl.plot(x[1],sim["BV_ms_2"],to_percentile(sim,"BV_hc_2","BV_ms_2"),
         marker='o',linewidth=0,label="MS 2",markersize=5,alpha=1,mew=0.5,color=ppl.almost_black)"""

ppl.scatter(x[1],sim["BV_hc_2"],to_percentile(sim,"BV_hc_2","BV_hc_2"),
         marker='o',label="HC 2",alpha=0.5, color=ppl.set2[0],s=200)
ppl.scatter(x[1],sim["BV_ms_2"],to_percentile(sim,"BV_hc_2","BV_ms_2"),
         marker='o',label="MS 2",alpha=1,color=ppl.almost_black,s=25)

Y2 = to_percentile(sim,"BV_hc_2","BV_euro_2")
#ppl.plot(x[1],sim["BV_euro_2"],Y2,
#         marker='o',linewidth=0,label="msHC 2",markersize=10,alpha=1,mew=0.5,color=ppl.set2[1])
scatter2 = ppl.scatter(x[1],sim["BV_euro_2"],Y2,
         marker='o',label="msHC 2",alpha=1,color=ppl.set2[1],s=100)

img_str = '%d: %2.2f, %2.2f'
labels = [img_str%(idx,X,Y2[idx]) for idx,X in enumerate(sim["BV_euro_2"])]
tooltip = plugins.PointLabelTooltip(scatter2, labels)
plugins.connect(fig, tooltip);

labels = [img_str%(idx,X,Y1[idx]) for idx,X in enumerate(sim["BV_euro_1"])]
tooltip = plugins.PointLabelTooltip(scatter1, labels)
plugins.connect(fig, tooltip);

x[0].set_ylim([-0.1,1.1])
x[1].set_ylim([-0.1,1.1])
x[0].set_ylabel("Percentile",size=14)
x[0].set_xlabel("Brain Volume",size=14)
x[0].set_title("Site 1",size=16)
x[1].set_title("Site 2",size=16)
x[1].set_ylabel("Percentile",size=14)
x[1].set_xlabel("Brain Volume",size=14)
#ppl.legend(x[0],loc=0,fontsize=14);
#ppl.legend(x[1],loc=0,fontsize=14);

In [66]:

pdiff = np.abs(to_percentile(sim,"BV_hc_2","BV_euro_2")-to_percentile(sim,"BV_hc_1","BV_euro_1"))
p,t,c = test_all(pdiff,linspace(0.01,0.25,40))

In [100]:

fig,x=subplots(1,3,figsize=(12,4))
m = np.mean(pdiff)
s = np.std(pdiff)
sns.kdeplot(pdiff,ax=x[0])
x[0].vlines(m,0,5)
x[0].vlines(m+s,0,5,linestyle='--',color="grey")
x[0].vlines(m-s,0,5,linestyle='--',color="grey")
x[0].set_xlim([0,0.3]);
x[0].set_ylim([0,5]);
x[0].set_title("Percentile Differences for MultiSite Controls",size=14)
x[1].set_title("Power Curve",size=14)
x[2].set_title("Equivalence Test Cutoff",size=14)
ppl.plot(x[2],linspace(0.01,0.25,40),t,color=ppl.set2[0],linewidth=3,linestyle="--",label="T")
ppl.plot(x[2],linspace(0.01,0.25,40),c,color=ppl.set2[1],linewidth=3,label="Cutoff")
ppl.legend(x[2],fontsize=14,loc=0)
powercurve = ppl.scatter(x[1],linspace(0.01,0.25,40)[:20],p[:20],marker='o',s=100);

img_str = '%2.2f, %2.2f'
labels = [img_str%(X,p[:20][idx]) for idx,X in enumerate(linspace(0.01,0.25,40)[:20])]
tooltip = plugins.PointLabelTooltip(powercurve, labels)
plugins.connect(fig, tooltip);

Aside: Running things in Parallel

In [56]:

from IPython.parallel.client.client import Client
from itertools import repeat
c = Client()
view = c[:]
view.block = True

In [27]:

def eq_test(y,ep):
    import rpy2.robjects as ro
    eq = ro.packages.importr("equivalence")
    import numpy as np
    x = ro.FloatVector(y)
    e = ep/np.std(y)
    res = eq.ptte_data(x,Epsilon=e)
    P = res.rx("Power")[0][0]
    T = res.rx("Tstat")[0][0]
    C = res.rx("cutoff")[0][0]
    return P,T,C

In [28]:

def test_all_parallel(pdiff,test_E = arange(0.001,0.01,0.0005),view=None):
    g = view.map(eq_test,[pdiff]*len(test_E),test_E)
    g = np.asarray(g)
    return g

Logistic Regression¶

We first run a logistic regression of the two sites independently: with brain volume and age and try to predict group membership:

Site 1:

In [29]:

import sklearn.linear_model as lm
a = lm.LogisticRegression()
X1 = np.vstack((sim["BV_ms_1"],sim["age_ms_1"])).T
y1 = sim["llr_1"]>0
X1 = sm.add_constant(X1)
a.fit(X1,y1)
acc1 = a.score(X1,y1)
print acc1

0.692

In [30]:

a = lm.LogisticRegression()
X2 = np.vstack((sim["BV_ms_2"],sim["age_ms_2"])).T
X2 = sm.add_constant(X2)
y2 = sim["llr_2"]>0
a.fit(X2,sim["llr_2"]>0)
acc2 = a.score(X2,sim["llr_2"]>0)
print acc2

0.676

In [64]:

_,x = subplots(1,2,figsize=(12,4))

sns.set_style("darkgrid")
df1 = pd.DataFrame(index=range(params["n_ms"][0]))
df1["pBV"] = sim["pBV_ms_1"]
df1["age"] = sim["age_ms_1"]
df1["EDSS"] = sim["dis_1"]
df1["Low Group"] = sim["llr_1"] > 0
iax = sns.interactplot("age","pBV","Low Group", df1,
                       cmap="coolwarm", filled=True, levels=25,ax=x[0],
                       colorbar=True,logistic=True,  
                       scatter_kws={"color": "black"},
                       contour_kws={"alpha": .8});
#iax.figure.colorbar(iax)

sns.set_style("darkgrid")
df2 = pd.DataFrame(index=range(params["n_ms"][1]))
df2["pBV"] = sim["pBV_ms_2"]
df2["age"] = sim["age_ms_2"]
df2["EDSS"] = sim["dis_2"]
df2["Low Group"] = sim["llr_2"] > 0
sns.interactplot("age","pBV","Low Group", df2,
                 cmap="coolwarm", filled=True, levels=25,ax=x[1],
                 colorbar=True,logistic=True,  
                 scatter_kws={"color": "black"},
                 contour_kws={"alpha": .8});
x[1].set_title("Site 2\nAccuracy = %2.2f"%(acc2),size=16);

x[0].set_title("Site 1\n Accuracy = %2.2f"%acc1,size=16);

We then calculate and plot the ROC curve for both sites:

In [33]:

def get_ROC(X,y):
    from sklearn.metrics import roc_curve, auc
    from sklearn.utils import shuffle
    import sklearn.linear_model as lm
    X, y = shuffle(X, y)
    half = int(len(y) / 2)
    X_train, X_test = X[:half], X[half:]
    y_train, y_test = y[:half], y[half:]

    # Run classifier
    classifier = lm.LogisticRegression()
    classifier.fit(X_train, y_train)
    probas_ = classifier.predict_proba(X_test)

    # Compute ROC curve and area the curve
    fpr, tpr, thresholds = roc_curve(y_test, probas_[:, 1])
    roc_auc = auc(fpr, tpr)
    return fpr,tpr,roc_auc

In [34]:

roc_info = view.map(get_ROC,[X1]*100,[y1]*100)
roc_info2 = view.map(get_ROC,[X2]*100,[y2]*100)
roc_info = np.asarray(roc_info)
roc_info2 = np.asarray(roc_info2)
fpr = roc_info[:,0]
tpr = roc_info[:,1]
auc = roc_info[:,2]
fpr2 = roc_info2[:,0]
tpr2 = roc_info2[:,1]
auc2 = roc_info2[:,2]

In [63]:

_,x=subplots(1,2,figsize=(12,4))
x[0].errorbar(np.mean(fpr),np.mean(tpr),xerr = np.std(fpr),yerr=np.std(tpr),alpha=0.15,color=ppl.set1[0],linewidth=3)
x[0].plot(np.mean(fpr),np.mean(tpr),color=ppl.set1[0],linewidth=3,label="Site 1")
x[0].errorbar(np.mean(fpr2),np.mean(tpr2),xerr = np.std(fpr2),yerr=np.std(tpr2),color=ppl.set1[1],alpha=0.15,linewidth=3)
x[0].plot(np.mean(fpr2),np.mean(tpr2),color=ppl.set1[1],linewidth=3,label="Site 2")
x[0].set_xlim([0,1])
x[0].set_ylim([0,1])
ppl.plot(x[0],[0, 1], [0, 1],color=ppl.set1[2],linewidth=2,linestyle='--')
x[1].set_title("Area Under the Curve",size=16)
x[0].set_title("ROC Curve",size=16)
x[0].set_xlabel("False Positive Rate")
x[0].set_xlabel("True Positive Rate")
x[0].legend(loc=0,fontsize=14)
ppl.boxplot(x[1],[auc,auc2]);
x[1].set_xticklabels(["Site 1","Site 2"]);

Next, we combine the data, run the regression and compute the ROC curve.

In [36]:

df_all = pd.DataFrame(index=range(sum(params["n_ms"])))
df_all["age"] = sim["age_all_ms"]
df_all["pBV"] = np.hstack((sim["pBV_ms_1"],sim["pBV_ms_2"]))
df_all["llr"] = np.hstack((sim["llr_1"],sim["llr_2"]))>0

In [37]:

a = lm.LogisticRegression()
X = np.vstack((df_all["pBV"],df_all["age"])).T
X = sm.add_constant(X)
y = df_all["llr"]
a.fit(X,y)
acc = a.score(X,y)

In [38]:

roc_info_all = view.map(get_ROC,[X]*100,[y]*100)
roc_info_all = np.asarray(roc_info_all)
fpr_all = roc_info_all[:,0]
tpr_all = roc_info_all[:,1]
auc_all = roc_info_all[:,2]

In [41]:

minlen = min([len(x) for x in fpr_all])
fpr_all = np.asarray([x[:minlen] for x in fpr_all])
tpr_all = np.asarray([x[:minlen] for x in tpr_all])

In [73]:

# Plot ROC curve
_,x=subplots(1,2,figsize=(12,4))

sns.interactplot("age","pBV","llr", df_all,cmap="coolwarm", 
                 filled=True, levels=20,colorbar=True,logistic=True,ax=x[0],
                 scatter_kws={"color": "black"},
                 contour_kws={"alpha": .8});
x[0].set_title("Logistic Regression\n %2.2f Accuracy"%(acc),size= 16)

x[1].errorbar(np.mean(fpr),np.mean(tpr),xerr = np.std(fpr),yerr=np.std(tpr),alpha=0.15,color=ppl.set1[0],linewidth=3)
x[1].plot(np.mean(fpr),np.mean(tpr),color=ppl.set1[0],linewidth=3,label="Site 1, AUC = %2.2f"%(np.mean(auc)))
x[1].errorbar(np.mean(fpr2),np.mean(tpr2),xerr = np.std(fpr2),yerr=np.std(tpr2),color=ppl.set1[1],alpha=0.15,linewidth=3)
x[1].plot(np.mean(fpr2),np.mean(tpr2),color=ppl.set1[1],linewidth=3,label="Site 2, AUC = %2.2f"%(np.mean(auc2)))

x[1].errorbar(np.mean(fpr_all,axis=0),np.mean(tpr_all,axis=0),xerr = np.std(fpr_all,axis=0),
              yerr=np.std(tpr_all,axis=0),alpha=0.15,color=ppl.set1[2],linewidth=3)
x[1].plot(np.mean(fpr_all,axis=0),np.mean(tpr_all,axis=0),color=ppl.set1[2],linewidth=3,
          label="Combined, AUC = %2.2f"%(np.mean(auc_all)))


ppl.plot(x[1],[0, 1], [0, 1],color=ppl.set1[3],linewidth=2,linestyle="--")
x[1].set_xlim([0.0, 1.0])
x[1].set_ylim([0.0, 1.0])
x[1].set_xlabel('False Positive Rate', size=14)
x[1].set_ylabel('True Positive Rate',size=14)
x[1].set_title('ROC',size=16)
x[1].legend(loc="lower right",fontsize=14);

Here we see that the ROC curve has improved with combined, calibrated data

Checking Calibration: Multisite Controls¶

We use the multisite controls to test the percentile calibration. We need to find what change in percentile brain volume to use as an equivalence margin \(\delta\) to distinguish biased HC populations.

In [43]:

params = {'n_hc':50, #number of controls
          'n_ms':[250,250], #number of MS patients at each site
          'age_hc_range':(20,52), # Age range of the controls
          'age_ms_range':[(20,45),(20,45)], # Age range of MS patients at each site
          'd_dist':[.65,.65], # Distribution weighting of EDSS scores for each site
          'alpha_a':-3, # Real coefficient of age
          'alpha_d':[-16,-16], # Real coefficient of EDSS score
          'b0':1644, #Brain Volume Intercept (in ccs)
          'sig_b':[57,57], # Brain variation across subjects (in ccs) for each site
          'sig_s':[16,16], # Variability in brain volume measurement across scanners for each site
          'alpha':[0.7,1.3], # Scaling of each scanner
          'gamma':[20,-20], # Offset of each scanner
          'lo_mean':[2,2], # EDSS distribution parameters at both sites
          'hi_mean':[6,6],  #EDSS distribution parameters at both sites
          'lo_sig': [1.1,1.1], #EDSS distribution parameters at both sites
          'hi_sig': [0.8,0.8], #EDSS distribution parameters at both sites
          }

In [53]:

b = boot(n_boot=100,keep_subs=False,params = params,age_euro_orig = age_euro,n_scans=1,sim_type=run_percentile_sim)
foo = [test_all_parallel(diff,linspace(0.01,0.2,20),view) for diff in b]
foo = np.asarray(foo)

In [55]:

_,x=subplots(1)
#pal = sns.blend_palette(["seagreen", "lightblue"])
sns.set_style("darkgrid")
pal = sns.color_palette("RdPu_r", 20)
sns.boxplot(foo[:,:,0],positions=range(20),ax=x,color=pal);
x.set_ylim([0,1.1])
x.set_title("Bootstrapped Power Analysis for Percentile Calibration",size=14)
x.set_xlabel("$\delta$ : Percentile",size=14)
x.set_ylabel("Power",size=14)
x.set_xticklabels(['%2.1f'%s for s in linspace(0.01,0.2,20)*100],rotation=90);

We only have the power to detect a 10% change in percentile across sites. This is really very high!

Effect of Scaling Parameter¶

In [44]:

a = deepcopy(params)
n_boot = 100
age_euro = np.array([51,32,24,30,32,37,39,24,24,37,29,32,40,41])
alpha_enum = [[0.9,1.1],[0.8,1.2],[0.7,1.3]]
accept_diff = linspace(0.05,0.2,30)
Pall = np.zeros((len(alpha_enum),n_boot,len(accept_diff)))
Tall = np.zeros((len(alpha_enum),n_boot,len(accept_diff)))
Call = np.zeros((len(alpha_enum),n_boot,len(accept_diff)))

for i,sig in enumerate(alpha_enum):
    a["alpha"] = sig
    b = boot(n_boot=n_boot,keep_subs=True,params = a,age_euro_orig = age_euro,n_scans=1,sim_type=run_percentile_sim)
    all_vals = np.asarray([test_all_parallel(diff,accept_diff,view) for diff in b])
    Pall[i,:,:] = all_vals[:,:,0]
    Tall[i,:,:] = all_vals[:,:,1]
    Call[i,:,:] = all_vals[:,:,2]

In [45]:

average_pass = np.mean(Call>Tall,axis=1)
vector = [x[1] - x[0] for x in alpha_enum]
_,x = subplots(1)
for i,v in enumerate(vector):
    ppl.plot(x,(np.array(accept_diff)+0.0*i)*100,average_pass[i],
             label="alpha difference = %2.2f"%(v),color=ppl.set2[i],marker='o',markersize=10);


x.legend(loc=0)
x.set_xlabel("Similarity Margin(Difference in Percentile Brain Volume)",size=14)
x.set_ylabel("Pass Rate/ Similarity Probability",size=14)
x.set_ylim([0,1.1])
x.set_title("Probability of Similarity ",size=16)
#x.set_xlim([1,5]);

Effect of Number of HC¶

In [46]:

a0 = deepcopy(params)
n_boot = 100
age_euro = np.array([51,32,24,30,32,37,39,24,24,37,29,32,40,41])
n_hc = [25,50,75,100, 150,500]
accept_diff =  linspace(0.05,0.2,30)
Pall = np.zeros((len(n_hc),n_boot,len(accept_diff)))
Tall = np.zeros((len(n_hc),n_boot,len(accept_diff)))
Call = np.zeros((len(n_hc),n_boot,len(accept_diff)))
t0 = time.time()
for i,hc in enumerate(n_hc):
    a0["n_hc"] = hc
    b = boot(n_boot=n_boot,keep_subs=True,params = a0,age_euro_orig = age_euro,n_scans=1,sim_type=run_percentile_sim)
    #b = boot(n_boot=n_boot,keep_subs=True,params = a0,age_euro_orig = age_euro,sim_type=run_percentile_sim,n_scans=1)
    all_vals = np.asarray([test_all_parallel(diff,accept_diff,view) for diff in b])
    Pall[i,:,:] = all_vals[:,:,0]
    Tall[i,:,:] = all_vals[:,:,1]
    Call[i,:,:] = all_vals[:,:,2]
    
t_total = time.time()-t0

In [47]:

average_pass = np.mean(Call>Tall,axis=1)
vector = n_hc
_,x = subplots(1)
for i,v in enumerate(vector):
    ppl.plot(x,(np.array(accept_diff)+0*i)*100,average_pass[i],
             label="n hc = %d"%(v),color=ppl.set2[i],marker='o',markersize=10);


x.legend(loc=0)
x.set_xlabel("Similarity Margin(Difference in Percentile Brain Volume)",size=14)
x.set_ylabel("Pass Rate/ Similarity Probability",size=14)
x.set_ylim([0,1.1])
x.set_title("Probability of Similarity ",size=16)
#x.set_xlim([0.1,6]);

Effect of Number of Scans¶

In [52]:

a = deepcopy(params)
n_boot = 100
age_euro = np.array([51,32,24,30,32,37,39,24,24,37,29,32,40,41])
n_scans = [1,2,3,4,5]
accept_diff = linspace(0.05,0.2,30)
Pall = np.zeros((len(n_scans),n_boot,len(accept_diff)))
Tall = np.zeros((len(n_scans),n_boot,len(accept_diff)))
Call = np.zeros((len(n_scans),n_boot,len(accept_diff)))
t0 = time.time()
for i,ns in enumerate(n_scans):
    b = boot(n_boot=n_boot,keep_subs=True,params = a,age_euro_orig = age_euro,n_scans=ns, sim_type=run_percentile_sim)
    all_vals = np.asarray([test_all_parallel(diff,accept_diff,view) for diff in b])
    Pall[i,:,:] = all_vals[:,:,0]
    Tall[i,:,:] = all_vals[:,:,1]
    Call[i,:,:] = all_vals[:,:,2]
    
t_total = time.time()-t0

In [53]:

average_pass = np.mean(Call>Tall,axis=1)
vector = n_scans
_,x = subplots(1)
for i,v in enumerate(vector):
    ppl.plot(x,(np.array(accept_diff)+0*i)*100,average_pass[i],
             label="n scans = %d"%(v),color=ppl.set2[i],marker='o',markersize=10);


x.legend(loc=0)
x.set_xlabel("Similarity Margin( Percentile Difference in Brain Volume)",size=14)
x.set_ylabel("Pass Rate/ Similarity Probability",size=14)
x.set_ylim([0,1.1])
x.set_title("Probability of Similarity ",size=16)
#x.set_xlim([0.1,6]);

In [57]:

Bias_All = []
age_euro = np.array([51,32,24,30,32,37,39,24,24,37,29,32,40,41])
a_hc = [[1,1],[1,1.01],[1,1.025],[1,1.04],[1,1.05]]
for idx,a in enumerate(a_hc):
    tmp = deepcopy(params)
    tmp["alpha_hc"] = a
    tmp["gamma_hc"] = [0, 0]
    b = boot(n_boot=150,keep_subs=False,params = tmp,age_euro_orig = age_euro,n_scans=2,sim_type=run_percentile_sim)
    foo1 = np.asarray([test_all_parallel(diff,linspace(0.05,0.4,50),view) for diff in b])
    Bias_All.append(foo1)

In [59]:

_,x=subplots(1,figsize=(8,6))
b = [0,1,2.5,4,5]
for idx, B in enumerate(Bias_All):
    ppl.plot(x,linspace(0.05,0.4,50)*100,np.mean(B[:,:,2]>B[:,:,1],axis=0),
             marker='o',label="%2.1f Percent Bias"%(b[idx]),color=ppl.set2[idx],markersize=10)


#x.set_xticks(range(20))
x.set_xlabel("Equivalence Margin $\delta$ (\% BV)",size=14)
x.set_ylabel("Probability to detect Similarity",size=14)
x.legend(loc=0)
x.set_ylim([0,1.1]);

In [60]:

B_0 = np.mean(Bias_All[0][:,:,2]>Bias_All[0][:,:,1],axis=0)
B_1 = np.mean(Bias_All[1][:,:,2]>Bias_All[1][:,:,1],axis=0)
B_2 = np.mean(Bias_All[2][:,:,2]>Bias_All[2][:,:,1],axis=0)
#B_3 = np.mean(Bias_All[3][:,:,2]>Bias_All[3][:,:,1],axis=0)
plot(B_1,B_0,label="1 percent bias, AUC = %2.2f"%np.trapz(B_0,B_1),color=ppl.set2[1],linewidth=4);
plot(B_2,B_0,label="2.5 percent bias, AUC = %2.2f"%np.trapz(B_0,B_2),color=ppl.set2[2],linewidth=4);
#plot(B_3,B_0,label="4 percent bias, AUC = %2.2f"%np.trapz(B_0,B_3),color=ppl.set2[3],linewidth=4);
plot([0,1],[0,1],linestyle='--',color=ppl.almost_black,alpha=0.5);
title("ROC Curve",size=16)
xlabel("False Positive Rate",size=14)
ylabel("True Positive Rate",size=14)
ylim([-0.05,1.05])
xlim([-0.05,1.05])
legend(loc=4,fontsize=14);

The percentile calibration method is more sensitive to bias than the OLS method, as shown by the area under the ROC curves for given levels of bias
With the percentile calibration, we still see that we need >50 HCs and >2 scans per multisite control
We have the power to detect a 10 percentile equivalence margin, but the data is only equivalent at margins > 15 percentiles

Percentile Calibration

Table of Contents

Background¶

Simulation¶

Assigning Groups: Log-Likelihood Ratio¶

\(LLR = log(\frac{Pr_{\theta_L}(E)}{Pr_{\theta_H}(E)})\)

Calibration by Percentile¶

Aside: Running things in Parallel

Logistic Regression¶

Checking Calibration: Multisite Controls¶

Power Curve¶

Effect of Scaling Parameter¶

Effect of Number of HC¶

Effect of Number of Scans¶

Effect of HC Bias¶

Conclusions¶

Next Steps:¶